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A THREE-DIMENSIONAL UNSTEADY CFD MODEL OF COMPRESSOR STABILITY 


Rodrick V. Chima 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 


ABSTRACT 

A three-dimensional unsteady CFD code called CSTALL has been 
developed and used to investigate compressor stability. The code 
solved the Euler equations through the entire annulus and all blade 
rows. Blade row turning, losses, and deviation were modeled using 
body force terms which required input data at stations between blade 
rows. The input data was calculated using a separate Navier-Stokes 
turbomachinery analysis code run at one operating point near stall, 
and was scaled to other operating points using overall characteristic 
maps. No information about the stalled characteristic was used. 
CSTALL was run in a 2-D throughflow mode for very fast calcula- 
tions of operating maps and estimation of stall points. Calculated pres- 
sure ratio characteristics for NASA stage 35 agreed well with 
experimental data, and results with inlet radial distortion showed the 
expected loss of range. CSTALL was also run in a 3-D mode to inves- 
tigate inlet circumferential distortion. Calculated operating maps for 
stage 35 with 120 degree distortion screens showed a loss in range and 
pressure rise. Unsteady calculations showed rotating stall with two 
part- span stall cells. The paper describes the body force formulation in 
detail, examines the computed results, and concludes with observa- 
tions about the code. 


NOMENCLATURE 

b blockage 

E, G, H inviscid fluxes 

e total energy 

F body force for turning 

/ body force for loss 

h, p, T, s enthalpy, pressure, temperature, entropy 
K Centrifugal, Coriolis, and blockage source terms 

m meridional directional 

m c corrected mass flow 

q vector of conserved variables 

t time 

V velocity 


x, 0, r cylindrical coordinates 

u, v, w cylindrical velocity components 

p, %, 8 relative flow, blade, and deviation angles 

A difference 

streamwise grid direction 
p density 

T throughflow time scale 

body force vector 

(]) body force scaling function 

ty turn turning function for deviation 

Q blade row angular velocity 

co loss coefficient 


Subscripts 

0 

1,2 

LE, TE 
t, s, 5, ref 
ex, °° 


stagnation state 

upstream, downstream 

leading edge, trailing edge 

turning, entropy, deviation, reference point 

grid exit, downstream of throttle 


Superscripts 

V relative velocity 

ss steady state 


INTRODUCTION 

Compressor stall and surge can have catastrophic consequences in 
aircraft, yet prediction of these phenomena remains as one of the 
major unsolved problems in turbomachinery. Many models of stall 
that give some insight into the phenomena have been developed but 
few models are capable of predicting stall onset, except perhaps for 
full computational fluid dynamic (CFD) analysis of the entire com- 
pressor which is beyond the scope of this paper. Existing models 
range from analytic models of rotating stall or compression system 
stability to 2-D or 3-D CLD models of compression systems. All mod- 
els require input of some information about compressor performance, 
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Figure 1 — Meridional grid used for NASA stage 35. 


usually equivalent to axisymmetric pressure rise and loss characteris- 
tics for the machine over its entire operating range, including stall. 
This information can be difficult to obtain accurately. This paper 
describes a 3-D unsteady CFD model of compressor stability that uses 
a new method for obtaining information about compressor perfor- 
mance directly from steady CFD calculations performed separately. 
First some existing models and their input requirements will be 
described. 

Pampreen [1] gives an excellent review of several analytic models 
of compressor stall. Early models predict the propagation velocity of 
stall cells. Later system stability models predict the behavior of a duct/ 
compressor/plenum/throttle system. These models usually assume 
that the axisymmetric static pressure rise characteristic of the com- 
pressor is known at all flow conditions, including stall. 

Takata and Nagano [2] developed a nonlinear analysis of rotating 
stall that modeled incompressible flow through a 2-D isolated blade 
row with infinitesimal pitch but finite chord. The static pressure loss 
and exit flow angle were assumed to be known functions of the inci- 
dence angle, and changes in loss were assumed to lag in time behind 
changes in incidence. The numerical results predicted stall inception, 
stall cell growth, and stall cell propagation velocity. Several authors 
have extended these ideas to 2-D compressible flows using modern 
CFD methods. Longley [3] presented an unsteady 2-D model of com- 
pressor stability and tested the model for a hypothetical four-stage 
compressor. Demargne and Longley [4] applied the model to four real 
compressors. Lindau and Owen [5] applied a similar model to NASA 
compressor stage 35. Both models used measured blade row charac- 
teristics extended over a large flow range as input, and used time-lag 
equations to update the body-force terms gradually. 

Inviscid axisymmetric throughflow models have long been used to 
model steady flows in turbomachinery. These models represent the 
axisymmetric average flow on a meridional surface. Since the equa- 
tions are averaged tangentially they cannot predict turning, and since 
they are inviscid they cannot predict losses. Body force terms are 
added to model these effects. Marble [6] derived a formulation for 
body forces required to produce given changes in swirl (rv 0 ) and 

entropy along a streamline. Stewart [7] described an iterative formula- 
tion for body forces required to turn the flow to a given angle, and 
used that formulation in an axisymmetric model of a full engine. Both 
of these formulations were used in the present work and are described 
later. 

By incorporating axisymmetric throughflow codes or ideas in 
unsteady CFD codes, several authors have developed codes for model- 
ing circumferential distortion and stall. Throughflow codes usually 


incorporate extensive correlations for blade row turning and loss, and 
thus provide rough predictions of blade row performance at all operat- 
ing points. The codes often require additional calibration for each new 
case. Escuret and Gamier [8] developed two unsteady codes for com- 
pressor instability. The first was a 2-D throughflow analysis that mod- 
eled stall and surge. The second was a 3-D stall model. Both codes 
were coupled to a SNECMA throughflow code to provide the blade 
row characteristics, and both used time-lag models to update the 
body-force terms. Hale, et al. [9] merged the HT0300 throughflow 
code and the 3-D NPARC code (now called WIND) to produce his 
TEACC code. That code was used to predict the effects of circumfer- 
ential distortion on steady performance of an isolated rotor. Hale, et al. 
[10] used the same code to predict the performance of a 3-stage mili- 
tary fan coupled to an F-16 forebody and inlet. 

Gong et al. [11] described a 3-D computational model for compres- 
sor instability. The model was applied to the GE low- speed (incom- 
pressible) research compressor and was used to simulate both modal 
and spike stall inception. Gong’s dissertation [12] extended the model 
to compressible flows and used it to show the steady response of 
NASA stage 35 to inlet distortion. In this work the body forces were 
modeled using two constant force coefficients similar to lift and drag 
coefficients. Hsiao, et al. [13] incorporated the compressible model in 
the WIND Navier-Stokes code and used the code to analyze a nacelle/ 
fan/exit guide vane/nozzle configuration. Here the body forces were 
interpolated directly from 3-D Navier-Stokes solutions of the rotor 
and guide vane onto the entire throughflow grid, so the model was 
limited to one operating point. 

In the present work a new CFD code called CSTALL was devel- 
oped to investigate compressor stability. The code was used to solve 
the unsteady 3-D Euler equations through a compressor using a finite- 
difference scheme. The computational grid shown in Fig. 1 resolved 
the blade planform but was evenly spaced in the circumferential direc- 
tion. The equations were treated as axisymmetric in blade passages 
(shown in color.) Marble’s formulation for turning and loss forces was 
used for steady calculations, and Stewart’s formulation for turning to a 
specified angle was used for unsteady calculations. Body force input 
data was calculated using a Navier-Stokes turbomachinery analysis 
code run at one operating point near stall, and was scaled to other 
operating points using normalized characteristic maps. No informa- 
tion about the stalled characteristic was used. The input data was easy 
to obtain and allowed the inviscid axisymmetric throughflow code to 
represent the blade rows with the same accuracy as the underlying 
Navier-Stokes calculations. CSTALL was used in a 2-D throughflow 
mode for fast calculations of operating maps and estimation of stall 
points, in a 3-D mode for calculations of circumferential distortion, 
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and in a 3-D unsteady mode to investigate rotating stall in a compres- 
sor stage. 

GOVERNING EQUATIONS 

The governing equations are the Euler equations written in cylindri- 
cal coordinates (x, 0, r) in a stationary frame of reference. 


dq dE ldG dH v , 
dt dx r00 
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The energy e is given by e = p^C v T+ i V j and the pressure p is 

given by p = (y - 1 )^e - V 2 j , where V 2 = u 2 + v 2 + w 2 . 

<f> is the total body force per unit mass and will be discussed in 
detail later. K is a source term that includes the centrifugal and Corio- 
lis forces, and derivatives of the blockage b , which is the fraction of 
the annulus open to the flow. To account for blade thickness, b can be 
i i N 

modeled as b - 0 - 0 I — . To account for endwall blockage, b can 
' s p '2n 

be varied by some function of axial distance. 

Axisymmetric Equations in Blade Passages 

The full 3-D Euler equations (1-3) are solved in the annular ducts 
outside of blade passages. Within blade passages there should be no 
communication with neighboring passages, but the uniform computa- 
tional grid used here does not resolve the solid blade shapes. Instead it 
is assumed that blades are infinitely thin and that the flow is axisym- 
metric on each meridional grid surface. The blockage term b 
described above is used to account for variation of the throughflow 
velocity due to finite blade thickness. 

Stator passages are modeled as axisymmetric in the absolute frame 
of reference and the 0-derivatives can be dropped immediately: 


1 dG A . . , 

- = 0 in stator passages. 


(4) 


Rotor passages are modeled as axisymmetric in the relative frame 
of reference. Equations (1-2) are transformed to the relative frame, 
all 0-derivatives are set to zero, then the equations are transformed 
back to the absolute frame. The result is that 


1 dG ~dq . 

^ -> in rotor passages, 

r o0 d0 


where Q is the blade row rotational speed. 


(5) 


BODY FORCES 

Ducts are modeled using the 3-D Euler equations with all body 
forces <f> = 0 . Blade passages are modeled using the axisymmetric 
Euler equations which give no information about turning or loss. That 
information is specified using the body force term <f> . The formula- 
tion used for <I> is based on the formulation developed by Marble [6] 
and outlined below. 

The body force array used in (1) may be written as 


<S> 


br\o, <J> x> $ e> <J> r , 


in blade passages. (6) 


The term r£2d> 0 used in the energy equation represents the work done 
by a rotor. 

The total force vector <I> can be written as <I> = F + f , where F is 

— > 

a turning (or lift) force normal to the relative velocity V ' , 


-> — > 

F • V = 0 

uF x + v'F 0 + wF r = 0 


-> — > 

and / is a loss (or drag) force parallel but opposed to V 


(7) 


f\\V' 


V*Ufr) = -l/l^ 


( 8 ) 


-> 

Y 


where V - v - r£l is the relative tangential velocity component. 

The components of these forces can be found by starting with the 
steady, axisymmetric Euler equations in non-conservative form. The 
equations can be written along a meridional streamline m using the 
following transformations: 


dm 2 - dx 2 + dr 2 

,,2 _ t /2 . , a ,2 


V m d m = ud x + wd r 


The momentum equations become simply 

x; P v m d m u + d xP = ®x 

0; p v m d m (rv) = r * 0 


(9) 


( 10 ) 


r: P V m d m W + d rP = 
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The 0 -momentum equation immediately gives the total tangential 
force, 


^0 - F e + /e - 


v md(rv) 

P r dm ’ 


(ID 


showing that the change of angular momentum along a streamline is 
proportional to the total applied tangential force. 

Loss forces are derived from the canonical equation of state, 
Tds = dh- dp / p . This equation is applied along a streamline and 
the pressure gradient is eliminated using the momentum equations, 
giving: 


V m d m h 0~ r£l % 


Tv m d m s 


+ u(F + f ) + v'(F 0 + fa) + w(F + f ) 


( 12 ) 


where h 0 is the total enthalpy. Equation (7) can be used to eliminate 
three terms on the right hand side, leaving 

v m d m h 0 ~ = Tv m d m s + ( u fx + v 'fQ + w fr) ( 13 ) 

In the absence of losses the right side of (13) equals zero since the 
entropy gradient d m s and the loss forces /• are both zero. The terms 
remaining on the left side of (13) give: 


Body Force Required to Turn to a Specified Angle 

It is often desirable to specify the relative flow angle P = % + 8 
within the blade row, where % is the blade angle and 5 is the devia- 
tion. There does not seem to be a way to specify the turning force <L 0 

explicitly in terms of P ; however, an iterative scheme similar to a 
scheme described by Stewart [7] can be used. The angles are shown in 
Fig. 2 and the relative velocity V' is shown misaligned with the 
desired flow angle P + % . A turning force is defined to be proportional 
to the difference between the desired and actual tangential velocity 
components, 



Figure 2 — Nomenclature used for deviation model. 


dh n 

Vm d^ = ra ^ = rQ ( F e + /e> ( 14 ) 

Equations (11) and (14) can be combined to give the differential form 
of the Euler turbine equation: 




(19) 


3*0 _ ^ d(rv) 

dm dm 


(15) 


which will be used later. 

Equation (14) says that the rate of increase of total enthalpy along a 
streamline equals the rate at which torque does work on the system. 
This is a general application of the first law of thermodynamics and is 
true whether losses are present or not. Thus it follows that the entropy 
variation along a streamline can found from the right side of (13): 


When the flow turns to the desired angle F 0 + 1 = F 0 and the iter- 
ation is converged. 

y turn * s a turn i n g function that varies the body force from zero at 

the leading edge to maximum at the trailing edge. Other turning func- 
tions based on the local blade angle were investigated, but they often 
resulted in non-physical supersonic throughflow velocities at flow 
rates below the experimental rotor choking flow. 


Tv m d m s = -(“/* + v ’/e + w / r ) = -f‘V' (16) 

Solving for the magnitude of the loss force \f\ , 


I/I 


„ V m ds 
~^d^l 


(17) 


At this point f x , / 0 , f r , and F 0 = <f> 0 - / 0 are all known. The 
meridional turning force F m is found using equation (7,) then 
F x and F r are found by assuming that F m is tangent to the stream- 
wise (§) grid lines which follow the hub and casing. 


F 


m 



J- 


2 2 
x/+ 


F = 




(18) 


Body Force Unknowns 

Either equation (1 1) or (19) can be used to find the magnitude of the 
total tangential turning force. For convenience equation (11) will be 
referred to as the turning model and equation (19) will be referred to 
as the deviation model. Equation (13) gives the magnitude of the loss 
force. 

The term d(rv)/dm in (11) cannot be calculated using the axisym- 
metric equations, and the terms ds/dm in (17) and 5 in (19) cannot 
be calculated from the Euler equations. Here these terms have been 
modeled using results from steady Navier-Stokes calculations done 
previously using the author’s CFD code SWIFT, see Chima [14, 15.] 
However, these terms could also be calculated from experimental data 
or design intent using the same formulation. The terms were approxi- 
mated as a combination of spanwise reference profiles evaluated 
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Figure 3 — Top: Spanwise profiles of A (rv) , A s/R , and deviation computed with SWIFT. 

Bottom: Maps of 1-D variation of A(rv) , A s/R , and deviation with corrected flow computed with SWIFT 


upstream and downstream of each blade row at a single reference flow For a stator the temperature ratio equals one, and it can be shown 

near stall, and a scaling function that varied with corrected flow. that A s/R « co , the total pressure loss coefficient. 


d(rv) 

dm 


A (rv) 
Am 


ty t (rh c ) 

ref 


ds 

dm 


As 

Am 



5 = 5 r e / + ( M" I c) 


( 20 ) 


In (20), the terms labeled refare the spanwise reference profiles. 
They were evaluated from a single SWIFT solution and stored for 
each blade at each spanwise location. The gradient terms (AA) were 
evaluated along streamwise grid lines, although technically they 
should be evaluated along streamlines. They were differenced 
between stations upstream and downstream of a blade. Since rv and 
s are convected along streamlines, the exact axial location of these 
stations is unimportant. Since the gradients were modeled as con- 
stants, the turning and loss vary linearly through the blade row. The 
deviation 8 re y was evaluated at blade trailing edges. 


For a rotor the term A(rv )/ Am can be calculated from Euler’s tur- 
bine equation (15.) For a stator it must be input directly. 

The entropy gradient was evaluated using the canonical equation of 
state written in terms of total conditions and integrated across a blade 
row: 


The input reference profiles are shown at the top of Fig. 3 (the 
SWIFT calculations are discussed later.) The turning profiles show 
that the rotor adds swirl to the flow with overturning at the endwalls, 
and that the stator removes most of the swirl. The entropy rise profiles 
show that the rotor is especially lossy at the tip, and that the cantile- 
vered stator is lossy at the hub. The deviation profiles are fairly uni- 
form except at the endwalls. Since these profiles were evaluated by 
averaging a 3-D Navier-Stokes solution, they inherently include 
losses, blockage, and over- or under-turning due to shock-boundary 
layer interaction, tip clearance flows, and secondary flows, all of 
which are important for modeling stall in tip-critical blades. 

In (20), the terms <\>(m c ) are normalized scaling functions used to 
scale the reference profiles to other flow rates. They are defined as: 


A (rv)(m c ) 

ty t (m) = = 

A (rv)(m c ) ref 

^s(^ c ) = [8(m c ) 


<M m c) 


A s(m c ) 

A S&cKef 


5 (m ) ,] 

c ref J midspan 


( 22 ) 


The overbar denotes a spanwise average, i.e., A(rv) ~ r 2 v 2 -r l v l 


As 

~R 



and As = s 2 -s 1 . The deviation scaling function was evaluated at 
(21) midspan since the spanwise average deviation is poorly defined. 
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The scaling functions vary with m c , the corrected mass flow at the 
blade leading edges. They are normalized such that they have no 
effect at the reference flow (m c ) • They were evaluated using 

SWIFT calculations of the stage speed line but could also be evaluated 
using experimental data or design intent. The functions were input as 
short tables for each blade row and interpolated linearly as needed. 
The functions were extrapolated linearly outside the range of the input 
data, but the mass flow was restricted to be less than the calculated 
choking value and greater than some stall cutoff. When the code was 
used to determine the stall point the cutoff was set very small. When it 
was used to calculate behavior in stall the cutoff was set to about 80 
percent of the computed stall flow. The input scaling functions are 
shown at the bottom of Fig. 3. The turning, entropy rise, and deviation 
all increase as flow rate decreases, as expected. Surprisingly, the 
entropy rise is linear with flow rate for both blade rows. 

For 3-D flows body forces were calculated independently on each 
meridional grid plane. References profiles were the same for all planes 
but scaling functions were determined for each plane based on the 
inlet corrected flow per unit area for that plane. In that way body 
forces could respond to circumferential variations in upstream condi- 
tions. Note that body forces are also functions of the local velocity and 
thermodynamic properties that are computed as part of the solution. 

Body Force Lag Model 

The body force terms described above are functions of the corrected 
mass flow calculated at blade leading edges and respond instanta- 
neously everywhere to changes in that flow. In reality the blade forces 
respond gradually to changes in flow or incidence. Time-lag models 
are commonly used to simulate this finite response rate, (see Longley 
[3] or Lindau, et al. [5].) The following time-lag model was used in 
the present work: 

= (F ss -F) (23) 

where t is a time scale associated with the blade response, typically 
the throughflow time, F is F 0 or |/| (the lag model is applied to 

both,) and F ss is the steady state value of the force. Computationally 
the model is applied using, 

F n+i = (l - r )F n +rF ss . (24) 

where r = At /% is an under relaxation factor. For steady calculations 
r — 1 . For the unsteady calculations shown later the throughflow 
time is about 100 A t , so that r = 0.01 . 

NUMERICAL MODEL 

The Euler equations (1-3) were transformed to body-fitted coordi- 
nates using finite-difference techniques described in Tweedt, et al. 
[16.] Since the computational grid was uniform in the tangential direc- 
tion, (A0 = constant) , the grid and metrics only required 2-D stor- 
age. The equations were discretized with either a central-difference 
scheme or the AUSM + upwind scheme described by Chima, et al. 
[15.] The equations were solved using an explicit Runge-Kutta 
scheme. For steady flows a two- stage scheme with spatially- varying 


time step and implicit residual smoothing was used to accelerate con- 
vergence. For unsteady flows a second-order four-stage scheme with 
constant time step was used. 

Boundary Conditions 

At the inlet the upstream-running Riemann invariant was extrapo- 
lated from the interior, and T 0 , and v were specified. In all cases 

T 0 was constant and v = 0 . The inlet P 0 profile reflected the turbu- 
lent inlet boundary layer profiles and resulting blockage used in 
SWIFT, but it was extrapolated to the wall on the coarser Euler grids 
used in CSTALL to avoid problems with zero velocity at the wall. For 
radial and circumferential distortion the inlet profile was modified 
locally as described later. 

At the exit, u, v, and w were extrapolated from the interior and p 
was set isentropically. The pressure was set using the radial equilib- 
rium equation and a throttle model, 



r hub 


Here K f is a dimensionless throttle coefficient discussed later and p ^ 
is the static pressure downstream of the throttle. 

RESULTS 

NASA Stage 35 

NASA stage 35 is a transonic inlet stage for a core compressor and 
was used here as a test case. The stage has a design pressure ratio of 
1.82 at a mass flow of 20.19 kg/sec and a rotor tip speed of 455 m/sec. 
The rotor has 36 multiple-circular- arc (MCA) blades with a maximum 
radius of 9.94 cm and a tip clearance of ~0.2 percent of the span. The 
cantilevered stator has 46 MCA blades with a hub clearance of -0.5 
percent span. The endwall clearances were much smaller than the grid 
spacing used in Fig. 1, but the effects of the clearances were captured 
by the reference profiles used in Fig. 3. Stage 35 was originally 
designed and tested at NASA Glenn Research Center by Reid and 
Moore [17.] Radial distributions of static and total pressure, total tem- 
perature, and flow angle were measured at axial stations located one 
rotor axial chord upstream of the rotor and 2.43 chords downstream of 
the stator. All calculations reported here were made at 100 percent 
design speed to allow comparison with this data. 

SWIFT - Navier-Stokes Analysis 

The SWIFT code is a multiblock Navier-Stokes analysis code for 
turbomachinery blade rows described by Chima, et al. [14, 15.] The 
code solves the Navier-Stokes equations on body-fitted grids using an 
explicit finite-difference scheme. It includes viscous terms in the 
blade-to-blade and hub-to-tip directions, but neglects them in the 
stream wise direction using the thin-layer approximation. The Bald- 
win-Lomax turbulence model was used for the present work. Multi 
blade row calculations were made using a steady mixing plane 
approach that transfers characteristic variables across the interface 
between blade rows. 
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A central difference scheme was used. Artificial dissipation and 
Eigenvalue scaling were used to capture shocks and to control point 
decoupling. An explicit, four-stage Runge-Kutta scheme was used to 
solve the flow equations. Calculations were run at a Courant number 
of 5.6 using a spatially-varying time step and implicit residual 
smoothing to accelerate convergence to a steady state. 

Calculations for stage 35 were run using a grid with five blocks: an 
H-block upstream, C-blocks around the rotor and stator, and O-blocks 
for the rotor tip clearance and stator hub clearance. Tip clearance grids 
had 9-11 points across the gap. The full grid had about 1.13 million 
points. The calculations were run several years ago on the Cray C-90 
at NASA Ames Research Center, on 6-8 processors. Cases near choke 
were run 3,000 iterations, which took about 1.1 hours wall clock time. 
Cases near stall were often run 3-4 times as long to ensure conver- 
gence. 

CSTALL - Euler Analysis 

CSTALL calculations were run on the axisymmetric grid shown in 
Fig. 1. The grid had 147 points axially and 21 points spanwise, for a 
total of 3,087 points. Steady 2-D cases were run with a 2-stage Runge- 
Kutta scheme on an SGI Octane2 with two processors. Cases near 
choke were converged to machine zero in 3,000 iterations, which took 
about 30 seconds. Cases near stall were often run 15-20,000 iterations 
to ensure convergence. 

3-D cases were run with 90 points evenly spaced in the circumfer- 
ential direction, for a total of 278,830 points. Steady cases near choke 
took about an hour. Unsteady 3-D cases were run with a 4-stage 
Runge-Kutta scheme with a constant time step and no residual 
smoothing. The time step was chosen to give a maximum Courant 
number of about 2.5. One rotor revolution required 2,500 time steps, 
which took about 40 minutes on the SGI. A dual time stepping 
scheme could be used to reduce CPU time for unsteady cases but was 
not used here. 
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Figure 4 — Experimental and computed characteristics of 
static-to-total pressure ratio for stage 35. 



Figure 5 — Spanwise distributions of experimental and 
computed total temperature ratio downstream of rotor 35. 


Predicted Performance 

Figure 4 compares static-to-total pressure ratio characteristics for 
stage 35 as measured by Reid and Moore [17] with characteristics 
computed by three different methods. The experimental data (solid 
black circles) show a non-zero slope near the stall point of 18.2 kg/ 
sec. The characteristic computed with the SWIFT code (plus signs) 
matches the experimental data closely but predicts stall at a higher 
flow of 19.23 kg/sec. The reason for the discrepancy in predicted stall 
point is unknown, but was found to be unaffected by changes in grid 
resolution, differencing scheme, and turbulence model. CSTALL 
solutions computed with the turning model (red squares) and the devi- 
ation model (black triangles) agree closely with each other and with 
the SWIFT solution up to the stall point. However, the CSTALL 
results predict stall at slightly lower flow than the experiment and a 
much lower flow than SWIFT. 

Since stall models typically require input of an axisymmetric char- 
acteristic that mimics stall, e.g. Gong [12,] it was not obvious in 
advance whether the CSTALL formulation could even predict stall. 
Subsequent results will show the mechanism responsible for stall in 
the calculations. For now, note that the CSTALL solutions were com- 
puted using the scaling functions computed with SWIFT (Fig. 3,) 
which were extrapolated linearly to the left of the predicted stall point 



Total Pressure Ratio 

Figure 6 — Spanwise distributions of experimental and 
computed total pressure ratio downstream of stator 35. 

as necessary. It was found that the predicted stall point depended 
strongly on the slopes of the scaling functions at the stall point. Small 
changes in the slopes could be used to make the predicted stall points 
match, but here the slopes were left as predicted by SWIFT. 
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In general the turning and deviation models gave similar results at 
stable operating points. However, at stalled operating points the turn- 
ing model gave non-physical swings in flow angles while the devia- 
tion model held the flow angle close to the blade angle even when the 
flow reversed direction. Hence the remainder of this paper will only 
show results from the turning model at stable points and from the 
deviation model at stalled points. 

CSTALL solutions were computed using reference profiles from 
the SWIFT solution at 19.74 kg/sec (labeled REF. on Fig. 4) and 
should agree with SWIFT at that flow rate. Spanwise profiles of total 
temperature ratio downstream of the rotor are compared in Fig. 5. The 
SWIFT predictions (black line) agree closely with experimental data 
(black circles) except near the tip, where a higher temperature ratio 
was predicted. The CSTALL predictions (red squares) agree reason- 
ably well with the SWIFT predictions, showing that the body forces 
for turning are behaving correctly. Spanwise profiles of total pressure 
ratio across the stator are compared in Fig. 6. The SWIFT predictions 
(black line) are about a point higher than the data (black circles,) but 
do show the large loss due to the stator hub clearance. The CSTALL 
predictions (red squares) agree reasonably well with the SWIFT pre- 
dictions, showing that the body forces for loss are also behaving cor- 
rectly. Minor differences between SWIFT and CSTALL in Ligs. 5 and 
6 could be due to a number of modelling simplifications used in 
CSTALL, such as estimation of the body forces using data along grid 
lines rather than streamlines, neglect of losses in ducts, etc. 

Stall Mechanism 

CSTALL calculations of a speed line (Fig. 4) behave like most 
Navier-Stokes calculations, with fast convergence near choke and 
slow convergence near stall. When the exit pressure is raised above 
the stall point the mass flow drops rapidly and large regions of reverse 
flow develop. Investigation of the flow field at the stall point suggests 
the mechanism by which the calculations predict stall. Figure 7 shows 
predicted axial velocity contours within the stage at the last stable 
operating point of Fig. 4. The plot shows a region of low axial velocity 
on the casing from the rotor trailing edge to the stator leading edge. 
This low-speed region is directly related to the high entropy rise 
across the rotor at the tip (Fig. 3, top center.) Figure 8 plots the axial 
velocity on the casing vs. axial distance. The upper plot is for the next 
to the last stable operating point on Fig. 4 and the lower plot is for the 
stall point. It is apparent that the axial velocity at the casing goes to 
zero quickly as the flow approaches stall, and that the stage stalls 
when the velocity reaches zero. 

Effect of Throttle Boundary Condition 

The static-to-total pressure ratio characteristics shown in Fig. 4 
were calculated using radial equilibrium at the exit (eq. 25) but not the 
throttle model, i.e., K f = 0. The throttle coefficient for the stage 35 

experiment was estimated to be K f ~ 1.314 using SWIFT results at 

the grid exit and assuming that the throttle exhausted back to the 
ambient. Another characteristic was calculated by holding 
K t = 1.314 constant and varying p ^ . The throttle characteristic at 

the stall point (black line) and the compressor characteristics pre- 
dicted with and without the throttle boundary condition are shown in 
Fig. 9. The throttle boundary condition reduced the predicted stall 



Figure 7 — Axial velocity component contours at stall. 



Axial position [cm] 

Figure 8 — Axial velocity along the casing of stage 35 at flow 
rates before stall (stable) and at stall. 
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Figure 9 — Effect of throttle exit boundary condition on the stall 
point. 


point from 17.1 kg/sec (red squares) to 16.4 kg/sec (blue triangles.) 
The slope of the static-to-total pressure ratio characteristics remained 
negative over the entire flow range. Since the CSTALL results predict 
a stall point slightly lower than the experiment even with K f = 0 , the 
throttle boundary condition was not used for subsequent calculations. 
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Inlet Total Pressure 

Figure 10 — Inlet total pressure profiles specified for stage 35. 



Figure 11 — Experimental and computed characteristics of 
static-to-total pressure ratio with and without radial distortion. 


Radial Distortion 

CSTALL requires input for a few operating points of a compressor. 
This input can be derived from steady Navier-Stokes calculations of 
periodic blade rows with clean inflow. Once this input is available 
CSTALL can be used to investigate the effects of inlet distortion on 
the performance of the compressor without the need for further 
Navier-Stokes computations. 

Radial distortion can have a big impact on compressor stability but 
it is easy to model with CSTALL. It was modeled by running 
CSTALL in a 2-D mode with the specified inlet total pressure profile 
modified as shown in Fig. 10. Figure 11 compares the predicted static- 
to-total pressure ratio characteristic with radial distortion to the com- 
puted and measured characteristics with clean flow shown previously. 
Radial distortion greatly reduces the actual mass flow and exit static 
pressure from their clean flow values, but the normalized characteris- 
tic of static-to-total pressure ratio versus corrected mass flow (red tri- 
angles) collapses onto the clean flow characteristic (blue squares.) The 
main effect of radial distortion is to increase the predicted stall flow 
from 17.1 kg/sec. to 18.8 kg/sec. 



Figure 12 — Experimental and computed characteristics of 
static-to-total pressure ratio for stage 35 with and without 
circumferential distortion. 



Figure 13 — Computed characteristic of static-to-total pressure 
ratio for stage 35 with circumferential distortion during stall and 
stall recovery. 

Circumferential Distortion 

Circumferential distortion is a 3-D problem that requires analysis of 
the entire annulus. This is costly to compute with a Navier-Stokes 
code but easy to compute with CSTALL because CSTALL uses an 
Euler formulation and a much coarser grid. Circumferential distortion 
was modeled by running CSTALL in a 3-D mode with the inlet total 
pressure profile reduced by 13 percent over a 120 degree segment of 
the inlet boundary. Figure 12 compares the predicted static-to-total 
pressure ratio characteristic with circumferential distortion to the 
computed and measured characteristics with clean flow shown previ- 
ously. Like radial distortion, circumferential distortion reduces the 
actual mass flow and exit static pressure from their clean flow values. 
Unlike radial distortion, the normalized static-to-total pressure ratio 
versus corrected mass flow characteristic does not collapse onto the 
clean characteristic, but becomes lower everywhere. 

At operating points near stall it was noted that many circumferential 
locations in the distorted region were well below the stall point on the 
body force calibration curves (Fig. 3.) In this case a parallel compres- 
sor model would predict that the stage would be stalled, but the 
present model predicts stable operation. 
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Figure 14 — Top: Inlet corrected mass flow for stage 35 during 
stall and stall recovery. 

Bottom: Casing static pressures at 12 circumferential locations 
at the leading edge of rotor 35 during stall and stall recovery. 

The calculations were pushed from the last stable operating point 
into stall by increasing the specified exit static pressure ratio from 
1.49 to 1.50. An unsteady rotating stall pattern developed gradually. 
After 16 rotor revolutions the exit pressure ratio was decreased to 1.10 
and the calculations were run another 14 revolutions to simulate stall 
recovery to a high flow point. Figure 13 shows the locus of the 
unsteady operating points on the static-to-total pressure ratio charac- 
teristics. During stall there was a large variation in inlet mass flow but 
little variation in pressure ratio which is effectively set by the inlet and 
exit boundary conditions. 

Figure 14 (top) shows the unsteady history of the inlet mass flow 
plotted against rotor revolutions. The increase in exit pressure reached 
the inlet boundary in about one revolution and the mass flow began to 
drop. An oscillation developed in the mass flow after several revolu- 
tions and became periodic after 10-12 revolutions. The exit pressure 
was decreased at 16 revolutions and the inlet mass flow began to 
increase about one revolution later. The oscillations are still evident at 
30 revolutions but are decaying asymptotically towards zero. The bot- 
tom of Fig. 14 shows unsteady casing pressure traces at the rotor lead- 
ing edge at 12 evenly spaced circumferential locations. Pressure 
oscillations grew gradually over about 10 revolutions then decayed 
slowly after the exit pressure was decreased at 16 revolutions. The 
four traces in the center are in the distorted flow region and show more 
variation in the signal. Analysis of the pressure traces shows that there 
are two stall sells that rotate at 0.94 Q . 

Figure 15 shows axial velocity contours on an unwrapped surface 
near the casing. Flow is from top to bottom, rotation is left to right, 
and the low total pressure distorted inflow is at the center. The upper 
figure shows a stable operating point where the distorted inflow is 
transported slightly in the direction of the rotation. The lower figure 
shows the flow during rotating stall where two stall cells are evident in 
the rotor. Since the tip speed is transonic these cells generate 


upstream-running shocks that produce the large pressure oscillations 
seen in Fig. 14. 

Figure 16 shows entropy contours on an annular surface near the 
rotor trailing edge. Rotation is counterclockwise. The left figure 
shows a stable operating point where the high entropy distorted inflow 
that started upstream at bottom center has been transported a few 
blade pitches counterclockwise. The right figure shows the flow dur- 
ing rotating stall. Animations have shown that the flow first becomes 
unsteady at the 4 o’clock position where the rotor leaves the distorted 
region. The flow remains quite unsteady at that location as two part- 
span stall cells form and rotate around the annulus. 

SUMMARY AND DISCUSSION 

A 3-D unsteady CFD code called CSTALL was developed and used 
to investigate distortion and stall phenomena in compressors. The 
code uses a body force formulation that can be calibrated easily from 
a few steady Navier-Stokes calculations. CSTALL was used in a 2-D 
throughflow mode to predict operating maps and estimate stall points 
quickly and to predict the effects of radial distortion. It was used in a 
3-D steady mode to predict the effects of circumferential distortion on 
the operating map. Finally, CSTALL was used in a 3-D unsteady 
mode to model rotating stall and recovery from stall. The following 
observations can be made about the code: 

• Marble’s formulation of body forces for turning and loss [6] pro- 
vided an elegant framework for deriving body force information 
from CFD solutions or experimental data. The body forces gave 
physically realistic spanwise distributions of turning, deviation, and 
loss, which are important for accurate prediction of tip-critical stall. 
Marble’s formulation was used for stable operating points, but an 
iterative scheme developed by Stewart [7] for turning the flow to a 
specified flow angle was found to be useful for stalled operating 
points. 

• Body force data was input at one operating point and was scaled to 
other operating points using normalized characteristic maps. This 
procedure gave good predictions of the entire operating line, but 
other formulations are possible. For example, data could input at 
many operating points and interpolated directly, or body force 
unknowns could be correlated with local incidence along the span. 
These formulations might improve predictions of stall but would 
complicate the input to the code. 

• CSTALL predicted compressor stall without specification of a 
stalled characteristic. The predicted stall point was at a slightly 
lower flow than the experiment and a much lower flow than Navier- 
Stokes calculations used for calibration. It was found that the pre- 
dicted stall point was sensitive to the slopes of the normalized char- 
acteristic maps used to scale the body force data to other mass 
flows. These maps were extrapolated linearly into stall, giving the 
correct trends for loss and deviation but probably not for A(rv) . 
Further calculations will be needed to find the best way to calculate 
the normalized characteristic maps and to extend them into stall. 

• A throttle exit boundary condition was shown to decrease the pre- 
dicted mass flow at which the compressor stalled. Since that mass 
flow was low even without the throttle condition, the throttle condi- 
tion was not used for most cases shown here. It would be interesting 
to see if this boundary condition can be used to stabilize Navier- 
Stokes calculations. 
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Figure 15 — Axial velocity contours on an unwrapped surface near the tip of stage 35. 
Top: Steady flow with circumferential distortion over 120 deg. 

Bottom: Stalled flow. 



Figure 16 — Entropy contours at the trailing edge of rotor 35. 
Left: Steady flow with circumferential distortion over 120 deg. 
Right: Stalled flow. 
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• One useful feature of CSTALL was its ability to predict the effects 
of inlet distortion on a compressor quickly without the need for fur- 
ther Navier-Stokes computations. Radial distortion was analyzed 
and was shown to reduce the stable flow range of the compressor 
without affecting the shape of the normalized operating characteris- 
tic. Circumferential distortion also reduced the stable flow range of 
the compressor, but it reduced the static-to-total pressure ratio at all 
flow rates as well. 

• Unsteady calculations showed the ability of the model to predict 
rotating stall inception and recovery. These calculations used a 
body force lag equation to model the delayed response of body 
forces to changes in flow. The model required a time lag parameter 
that is usually taken to be the blade row throughflow time. Unsteady 
calculations seemed to be insensitive to changes in this parameter, 
but further calculations are needed to determine its full effect. It 
may also be possible to calibrate this parameter using 2-D unsteady 
Navier-Stokes cascade calculations. 

• Further calculations are planned to validate CSTALL predictions of 
rotating stall and to investigate casing injection flow control for stall 
prevention and recovery. 
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